Antoine Godin (Kingston University)
Agent-based and stock-flow consistent modelling: theory and applications - Paris - July 17
You need to install all the required libraries This is for traditional libraries
install.packages("expm")
install.packages("igraph")For non-conventional libraries, such as the one need to visualize Direct Acyclical Graphs (DAG), you need to do the following
source("http://bioconductor.org/biocLite.R")
biocLite("Rgraphviz")Finally you can then download the PKSFC package from github and install it locally
install.packages("path/PKSFC_1.5.tar.gz", repos = NULL,
type="source")Now we’re ready to load the package:
library(PKSFC)## Loading required package: expm
## Loading required package: Matrix
##
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
##
## expm
## Loading required package: igraph
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
sim<-sfc.model("Models/SIM.txt",modelName="SIMplest model
from chapter 3 of Godley and Lavoie (2007)")datasim<-simulate(sim)matplot(sim$time,datasim$baseline[,c("Yd","C")],type="l",
xlab="", ylab="",lwd=2,lty=c(2,3))
lines(sim$time,vector(length=length(sim$time))
+datasim$baseline["2010","C"], lwd=2)
legend(x=1970,y=50,lwd=2, legend=c("Disposable Income",
"Consumption","Steady State"),lty=c(2,3,1), bty="n")#1. EQUATIONS
C = C
G = G
T = T
N = N
Yd = W*N - T
T = theta*W*N
C = alpha1*Yd + alpha2*H_h(-1)
H = H(-1) + G - T
H_h = H_h(-1) + Yd - C
Y = C + G
N = Y/W
#2. PARAMETERS
alpha1=0.6
alpha2=0.4
theta=0.2
#EXOGENOUS
G=20
W=1
#INITIAL VALUES
H=0
H_h=0
#3. Timeline
timeline 1945 2010#)(100>10) will return 1.(-x) where x is the lag.# character at the beginning of the line. Each comment exactly above an equation will be considered as the description of the equation and will be stored in the internal representation of the sfc model object.print(sim)## $name
## [1] "SIMplest model \n\t\tfrom chapter 3 of Godley and Lavoie (2007)"
##
## $simulated
## [1] FALSE
##
## $variables
## name initial value description
## [1,] "Yd" NA "1. EQUATIONS"
## [2,] "T" NA ""
## [3,] "C" NA ""
## [4,] "H" "0" ""
## [5,] "H_h" "0" ""
## [6,] "Y" NA ""
## [7,] "N" NA ""
## [8,] "alpha1" "0.6" "2. PARAMETERS"
## [9,] "alpha2" "0.4" ""
## [10,] "theta" "0.2" ""
## [11,] "G" "20" "EXOGENOUS"
## [12,] "W" "1" ""
##
## $endogenous
## name initial value lag description
## [1,] "Yd" NA "0" "1. EQUATIONS"
## [2,] "T" NA "0" ""
## [3,] "C" NA "0" ""
## [4,] "H" "0" "1" ""
## [5,] "H_h" "0" "1" ""
## [6,] "Y" NA "0" ""
## [7,] "N" NA "0" ""
##
## $equations
## endogenous value equation description
## [1,] "Yd" "W*N-T" "1. EQUATIONS"
## [2,] "T" "theta*W*N" ""
## [3,] "C" "alpha1*Yd+alpha2*H_h_1" ""
## [4,] "H" "H_1+G-T" ""
## [5,] "H_h" "H_h_1+Yd-C" ""
## [6,] "Y" "C+G" ""
## [7,] "N" "Y/W" ""
##
## $time
## [1] 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958
## [15] 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972
## [29] 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986
## [43] 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000
## [57] 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010
##
## $matrix
## Yd T C H H_h Y N
## Yd 0 1 0 0 0 0 1
## T 0 0 0 0 0 0 1
## C 1 0 0 0 0 0 0
## H 0 1 0 0 0 0 0
## H_h 1 0 1 0 0 0 0
## Y 0 0 1 0 0 0 0
## N 0 0 0 0 0 1 0
##
## $blocks
## $blocks[[1]]
## [1] 2 3 4 6 7 1 5
##
##
## attr(,"class")
## [1] "sfc"
scenario_i
| Y | d | T | C | H H_ | h | Y | N alpha1 | alpha2 | theta | G W | it | er bl | ock 1 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1945 | NA | NA | NA | 0.00000 | 0.00000 | NA | NA | 0.6 | 0.4 | 0.2 | 20 | 1 | 0 |
| 1946 | 30.76923 | 7.692308 | 18.46154 | 12.30769 | 12.30769 | 38.46154 | 38.46154 | 0.6 | 0.4 | 0.2 | 20 | 1 | 293 |
| 1947 | 38.34320 | 9.585799 | 27.92899 | 22.72189 | 22.72189 | 47.92899 | 47.92899 | 0.6 | 0.4 | 0.2 | 20 | 1 | 269 |
| 1948 | 44.75193 | 11.187984 | 35.93992 | 31.53391 | 31.53391 | 55.93992 | 55.93992 | 0.6 | 0.4 | 0.2 | 20 | 1 | 263 |
| 1949 | 50.17471 | 12.543678 | 42.71839 | 38.99023 | 38.99023 | 62.71839 | 62.71839 | 0.6 | 0.4 | 0.2 | 20 | 1 | 258 |
| 1950 | 54.76322 | 13.690805 | 48.45402 | 45.29943 | 45.29943 | 68.45402 | 68.45402 | 0.6 | 0.4 | 0.2 | 20 | 1 | 254 |
| 1951 | 58.64580 | 14.661450 | 53.30725 | 50.63798 | 50.63798 | 73.30725 | 73.30725 | 0.6 | 0.4 | 0.2 | 20 | 1 | 251 |
| 1952 | 61.93106 | 15.482766 | 57.41383 | 55.15521 | 55.15521 | 77.41383 | 77.41383 | 0.6 | 0.4 | 0.2 | 20 | 1 | 247 |
| 1953 | 64.71090 | 16.177725 | 60.88862 | 58.97749 | 58.97749 | 80.88862 | 80.88862 | 0.6 | 0.4 | 0.2 | 20 | 1 | 244 |
| 1954 | 67.06307 | 16.765767 | 63.82884 | 62.21172 | 62.21172 | 83.82884 | 83.82884 | 0.6 | 0.4 | 0.2 | 20 | 1 | 242 |
| 1955 | 69.05337 | 17.263341 | 66.31671 | 64.94838 | 64.94838 | 86.31671 | 86.31671 | 0.6 | 0.4 | 0.2 | 20 | 1 | 239 |
| 1956 | 70.73746 | 17.684366 | 68.42183 | 67.26401 | 67.26401 | 88.42183 | 88.42183 | 0.6 | 0.4 | 0.2 | 20 | 1 | 236 |
| 1957 | 72.16247 | 18.040617 | 70.20309 | 69.22339 | 69.22339 | 90.20309 | 90.20309 | 0.6 | 0.4 | 0.2 | 20 | 1 | 234 |
| 1958 | 73.36824 | 18.342061 | 71.71030 | 70.88133 | 70.88133 | 91.71030 | 91.71030 | 0.6 | 0.4 | 0.2 | 20 | 1 | 231 |
| 1959 | 74.38851 | 18.597128 | 72.98564 | 72.28421 | 72.28421 | 92.98564 | 92.98564 | 0.6 | 0.4 | 0.2 | 20 | 1 | 229 |
| 1960 | 75.25182 | 18.812955 | 74.06477 | 73.47125 | 73.47125 | 94.06477 | 94.06477 | 0.6 | 0.4 | 0.2 | 20 | 1 | 227 |
| 1961 | 75.98231 | 18.995577 | 74.97789 | 74.47567 | 74.47567 | 94.97789 | 94.97789 | 0.6 | 0.4 | 0.2 | 20 | 1 | 224 |
| 1962 | 76.60041 | 19.150104 | 75.75052 | 75.32557 | 75.32557 | 95.75052 | 95.75052 | 0.6 | 0.4 | 0.2 | 20 | 1 | 222 |
| 1963 | 77.12343 | 19.280857 | 76.40428 | 76.04471 | 76.04471 | 96.40428 | 96.40428 | 0.6 | 0.4 | 0.2 | 20 | 1 | 220 |
| 1964 | 77.56598 | 19.391494 | 76.95747 | 76.65322 | 76.65322 | 96.95747 | 96.95747 | 0.6 | 0.4 | 0.2 | 20 | 1 | 218 |
| 1965 | 77.94044 | 19.485111 | 77.42555 | 77.16811 | 77.16811 | 97.42555 | 97.42555 | 0.6 | 0.4 | 0.2 | 20 | 1 | 215 |
| 1966 | 78.25730 | 19.564324 | 77.82162 | 77.60378 | 77.60378 | 97.82162 | 97.82162 | 0.6 | 0.4 | 0.2 | 20 | 1 | 213 |
| 1967 | 78.52541 | 19.631351 | 78.15676 | 77.97243 | 77.97243 | 98.15676 | 98.15676 | 0.6 | 0.4 | 0.2 | 20 | 1 | 211 |
| 1968 | 78.75227 | 19.688067 | 78.44033 | 78.28437 | 78.28437 | 98.44033 | 98.44033 | 0.6 | 0.4 | 0.2 | 20 | 1 | 209 |
| 1969 | 78.94423 | 19.736056 | 78.68028 | 78.54831 | 78.54831 | 98.68028 | 98.68028 | 0.6 | 0.4 | 0.2 | 20 | 1 | 207 |
| 1970 | 79.10665 | 19.776663 | 78.88332 | 78.77165 | 78.77165 | 98.88332 | 98.88332 | 0.6 | 0.4 | 0.2 | 20 | 1 | 205 |
| 1971 | 79.24409 | 19.811023 | 79.05511 | 78.96062 | 78.96062 | 99.05511 | 99.05511 | 0.6 | 0.4 | 0.2 | 20 | 1 | 202 |
| 1972 | 79.36038 | 19.840096 | 79.20048 | 79.12053 | 79.12053 | 99.20048 | 99.20048 | 0.6 | 0.4 | 0.2 | 20 | 1 | 200 |
| 1973 | 79.45879 | 19.864697 | 79.32348 | 79.25583 | 79.25583 | 99.32348 | 99.32348 | 0.6 | 0.4 | 0.2 | 20 | 1 | 198 |
| 1974 | 79.54205 | 19.885513 | 79.42756 | 79.37032 | 79.37032 | 99.42756 | 99.42756 | 0.6 | 0.4 | 0.2 | 20 | 1 | 196 |
| 1975 | 79.61250 | 19.903126 | 79.51563 | 79.46719 | 79.46719 | 99.51563 | 99.51563 | 0.6 | 0.4 | 0.2 | 20 | 1 | 194 |
| 1976 | 79.67212 | 19.918030 | 79.59015 | 79.54916 | 79.54916 | 99.59015 | 99.59015 | 0.6 | 0.4 | 0.2 | 20 | 1 | 192 |
| 1977 | 79.72256 | 19.930640 | 79.65320 | 79.61852 | 79.61852 | 99.65320 | 99.65320 | 0.6 | 0.4 | 0.2 | 20 | 1 | 190 |
| 1978 | 79.76524 | 19.941311 | 79.70656 | 79.67721 | 79.67721 | 99.70656 | 99.70656 | 0.6 | 0.4 | 0.2 | 20 | 1 | 187 |
| 1979 | 79.80136 | 19.950340 | 79.75170 | 79.72687 | 79.72687 | 99.75170 | 99.75170 | 0.6 | 0.4 | 0.2 | 20 | 1 | 185 |
| 1980 | 79.83192 | 19.957980 | 79.78990 | 79.76889 | 79.76889 | 99.78990 | 99.78990 | 0.6 | 0.4 | 0.2 | 20 | 1 | 183 |
| 1981 | 79.85778 | 19.964445 | 79.82222 | 79.80445 | 79.80445 | 99.82222 | 99.82222 | 0.6 | 0.4 | 0.2 | 20 | 1 | 181 |
| 1982 | 79.87966 | 19.969915 | 79.84957 | 79.83453 | 79.83453 | 99.84957 | 99.84957 | 0.6 | 0.4 | 0.2 | 20 | 1 | 179 |
| 1983 | 79.89817 | 19.974543 | 79.87272 | 79.85999 | 79.85999 | 99.87272 | 99.87272 | 0.6 | 0.4 | 0.2 | 20 | 1 | 177 |
| 1984 | 79.91384 | 19.978460 | 79.89230 | 79.88153 | 79.88153 | 99.89230 | 99.89230 | 0.6 | 0.4 | 0.2 | 20 | 1 | 175 |
| 1985 | 79.92709 | 19.981774 | 79.90887 | 79.89975 | 79.89975 | 99.90887 | 99.90887 | 0.6 | 0.4 | 0.2 | 20 | 1 | 173 |
| 1986 | 79.93831 | 19.984578 | 79.92289 | 79.91518 | 79.91518 | 99.92289 | 99.92289 | 0.6 | 0.4 | 0.2 | 20 | 1 | 170 |
| 1987 | 79.94780 | 19.986950 | 79.93475 | 79.92823 | 79.92823 | 99.93475 | 99.93475 | 0.6 | 0.4 | 0.2 | 20 | 1 | 168 |
| 1988 | 79.95583 | 19.988958 | 79.94479 | 79.93927 | 79.93927 | 99.94479 | 99.94479 | 0.6 | 0.4 | 0.2 | 20 | 1 | 166 |
| 1989 | 79.96263 | 19.990657 | 79.95328 | 79.94861 | 79.94861 | 99.95328 | 99.95328 | 0.6 | 0.4 | 0.2 | 20 | 1 | 164 |
| 1990 | 79.96838 | 19.992094 | 79.96047 | 79.95652 | 79.95652 | 99.96047 | 99.96047 | 0.6 | 0.4 | 0.2 | 20 | 1 | 162 |
| 1991 | 79.97324 | 19.993310 | 79.96655 | 79.96321 | 79.96321 | 99.96655 | 99.96655 | 0.6 | 0.4 | 0.2 | 20 | 1 | 160 |
| 1992 | 79.97736 | 19.994340 | 79.97170 | 79.96887 | 79.96887 | 99.97170 | 99.97170 | 0.6 | 0.4 | 0.2 | 20 | 1 | 158 |
| 1993 | 79.98084 | 19.995210 | 79.97605 | 79.97366 | 79.97366 | 99.97605 | 99.97605 | 0.6 | 0.4 | 0.2 | 20 | 1 | 156 |
| 1994 | 79.98379 | 19.995947 | 79.97974 | 79.97771 | 79.97771 | 99.97974 | 99.97974 | 0.6 | 0.4 | 0.2 | 20 | 1 | 153 |
| 1995 | 79.98628 | 19.996571 | 79.98285 | 79.98114 | 79.98114 | 99.98285 | 99.98285 | 0.6 | 0.4 | 0.2 | 20 | 1 | 151 |
| 1996 | 79.98839 | 19.997098 | 79.98549 | 79.98404 | 79.98404 | 99.98549 | 99.98549 | 0.6 | 0.4 | 0.2 | 20 | 1 | 149 |
| 1997 | 79.99018 | 19.997545 | 79.98772 | 79.98650 | 79.98650 | 99.98772 | 99.98772 | 0.6 | 0.4 | 0.2 | 20 | 1 | 147 |
| 1998 | 79.99169 | 19.997923 | 79.98961 | 79.98857 | 79.98857 | 99.98961 | 99.98961 | 0.6 | 0.4 | 0.2 | 20 | 1 | 145 |
| 1999 | 79.99297 | 19.998242 | 79.99121 | 79.99033 | 79.99033 | 99.99121 | 99.99121 | 0.6 | 0.4 | 0.2 | 20 | 1 | 143 |
| 2000 | 79.99405 | 19.998513 | 79.99256 | 79.99182 | 79.99182 | 99.99256 | 99.99256 | 0.6 | 0.4 | 0.2 | 20 | 1 | 141 |
| 2001 | 79.99497 | 19.998741 | 79.99371 | 79.99308 | 79.99308 | 99.99371 | 99.99371 | 0.6 | 0.4 | 0.2 | 20 | 1 | 139 |
| 2002 | 79.99574 | 19.998935 | 79.99468 | 79.99414 | 79.99414 | 99.99468 | 99.99468 | 0.6 | 0.4 | 0.2 | 20 | 1 | 136 |
| 2003 | 79.99640 | 19.999099 | 79.99549 | 79.99504 | 79.99504 | 99.99549 | 99.99549 | 0.6 | 0.4 | 0.2 | 20 | 1 | 134 |
| 2004 | 79.99695 | 19.999237 | 79.99619 | 79.99581 | 79.99581 | 99.99619 | 99.99619 | 0.6 | 0.4 | 0.2 | 20 | 1 | 132 |
| 2005 | 79.99742 | 19.999355 | 79.99677 | 79.99645 | 79.99645 | 99.99677 | 99.99677 | 0.6 | 0.4 | 0.2 | 20 | 1 | 130 |
| 2006 | 79.99782 | 19.999454 | 79.99727 | 79.99700 | 79.99700 | 99.99727 | 99.99727 | 0.6 | 0.4 | 0.2 | 20 | 1 | 128 |
| 2007 | 79.99815 | 19.999538 | 79.99769 | 79.99746 | 79.99746 | 99.99769 | 99.99769 | 0.6 | 0.4 | 0.2 | 20 | 1 | 126 |
| 2008 | 79.99844 | 19.999609 | 79.99805 | 79.99785 | 79.99785 | 99.99805 | 99.99805 | 0.6 | 0.4 | 0.2 | 20 | 1 | 124 |
| 2009 | 79.99868 | 19.999669 | 79.99835 | 79.99818 | 79.99818 | 99.99835 | 99.99835 | 0.6 | 0.4 | 0.2 | 20 | 1 | 122 |
| 2010 | 79.99888 | 19.999720 | 79.99860 | 79.99846 | 79.99846 | 99.99860 | 99.99860 | 0.6 | 0.4 | 0.2 | 20 | 1 | 119 |
Principle: Solving \(Ax=b,\, A\in\mathbb{R}^{n\times n},\, b\in \mathbb{R}^n\) via an iterative algorithm, where each iteration can be represented by \(L x^{k+1} = b-Ux^{k},\, A=L+U\). Where \(L\) is lower triangular and \(U\) is upper triangular.
See Fenell et. al (2016)
simex<-sfc.model("Models/SIMEX.txt",modelName="SIMplest model with expectation")
layout(matrix(c(1,2),1,2))
plot.dag(sim,main="SIM" )
plot.dag(simex,main="SIMEX" )ch6 <- sfc.model("Models/ch6.txt",modelName="Chapter6_openmodel")
graphs = generate.DAG.collapse( adjacency = ch6$matrix )
par(mfrow = c(1,3))
# first plot the orgianl grpah
plot_graph_hierarchy( graphs$orginal_graph, main = "orginal graph" )
# plot hte nodes that form the strongly connected compoent
plot_graph_hierarchy( graphs$SCC_graph, main = "SCC nodes" )
# plot the result ing DAG when we take the condensation of the graph
plot_graph_hierarchy( graphs$DAG, main = "resulting DAG" )#Doing the simulations
datasimex<-simulate(simex)
init = datasimex$baseline[66,]
simex<-sfc.addScenario(simex,"G",25,1960,2010,init)
datasimex<-simulate(simex)
datasim<-simulate(sim)
init = datasim$baseline[66,]
sim<-sfc.addScenario(sim,"G",25,1960,2010,init)
datasim<-simulate(sim)
#Plotting it all
layout(matrix(c(1,2),1,2))
matplot(sim$time,datasim$scenario_1[,c("H","C","Yd")],
lty=1:3,type="l",xlab="",ylab="",main="SIM")
legend(x=1944,y=130,legend=c("Wealth","Consumption",
"Disposable Income"), lty=1:3,bty="n")
matplot(simex$time,datasimex$scenario_1[,c("H","C","Yd",
"Yd_e")],lty=1:4,type="l",xlab="",ylab="",main="SIMEX")
legend(x=1944,y=130,legend=c("Wealth","Consumption","Disposable Income",
"Expecetd Disposable Income"),lty=1:4,bty="n")Let’s see how much time it takes to run sim:
ptm <- proc.time()
data1<-simulate(sim)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],
"seconds"))## [1] "Elapsed time is 3.67 seconds"
Now lets play with some of the parameters of the simulate function:
ptm <- proc.time()
data2<-simulate(sim,tolValue = 1e-3)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],"seconds"))## [1] "Elapsed time is 0.176000000000002 seconds"
ptm <- proc.time()
data3<-simulate(sim, maxIter=10)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],"seconds"))## [1] "Elapsed time is 0.173999999999999 seconds"
kable(round(t(cbind(data1$baseline[c(1,2,20,40,66),c("Y")],
data2$baseline[c(1,2,20,40,66),c("Y")],
data3$baseline[c(1,2,20,40,66),c("Y")])),digits=3))| 1945 | 1946 | 1964 | 1984 | 2010 |
|---|---|---|---|---|
| NA | 38.462 | 96.957 | 99.892 | 99.999 |
| NA | 38.464 | 96.874 | 99.577 | 99.911 |
| NA | 34.006 | 94.965 | 99.672 | 99.991 |
The order of equations matters, if first compute variables that do not depend on current period, this speeds the process. Define blocks of equation independent from the others.
print(simex$blocks)## [[1]]
## [1] 7
##
## [[2]]
## [1] 3
##
## [[3]]
## [1] 6
##
## [[4]]
## [1] 8
##
## [[5]]
## [1] 2
##
## [[6]]
## [1] 1 4
##
## [[7]]
## [1] 5
ptm <- proc.time()
dataex<-simulate(simex,tolValue = 1e-10)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3]
,"seconds"))## [1] "Elapsed time is 0.117999999999999 seconds"
kable(round(t(dataex$baseline[c(1,2,20,40,66),c("G","Y",
"T","Yd","Yd_e","C","H","H_h")]),digits=3))| 1945 | 1946 | 1964 | 1984 | 2010 | |
|---|---|---|---|---|---|
| G | 20 | 20 | 20.000 | 20.000 | 20 |
| Y | NA | 20 | 98.559 | 99.983 | 100 |
| T | NA | 4 | 19.712 | 19.997 | 20 |
| Yd | 0 | 16 | 78.847 | 79.987 | 80 |
| Yd_e | NA | 0 | 78.559 | 79.983 | 80 |
| C | NA | 0 | 78.559 | 79.983 | 80 |
| H | 0 | 16 | 78.847 | 79.987 | 80 |
| H_h | 0 | 16 | 78.847 | 79.987 | 80 |
kable(dataex$baseline)| Yd | T | C | H | H_h | Y | Yd_e | N | alpha1 | alpha2 | theta | G | W | iter block 1 | iter block 2 | iter block 3 | iter block 4 | iter block 5 | iter block 6 | iter block 7 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1945 | 0.00000 | NA | NA | 0.00000 | 0.00000 | NA | NA | NA | 0.6 | 0.4 | 0.2 | 20 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1946 | 16.00000 | 4.00000 | 0.00000 | 16.00000 | 16.00000 | 20.00000 | 0.00000 | 20.00000 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1947 | 28.80000 | 7.20000 | 16.00000 | 28.80000 | 28.80000 | 36.00000 | 16.00000 | 36.00000 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1948 | 39.04000 | 9.76000 | 28.80000 | 39.04000 | 39.04000 | 48.80000 | 28.80000 | 48.80000 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1949 | 47.23200 | 11.80800 | 39.04000 | 47.23200 | 47.23200 | 59.04000 | 39.04000 | 59.04000 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1950 | 53.78560 | 13.44640 | 47.23200 | 53.78560 | 53.78560 | 67.23200 | 47.23200 | 67.23200 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1951 | 59.02848 | 14.75712 | 53.78560 | 59.02848 | 59.02848 | 73.78560 | 53.78560 | 73.78560 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1952 | 63.22278 | 15.80570 | 59.02848 | 63.22278 | 63.22278 | 79.02848 | 59.02848 | 79.02848 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1953 | 66.57823 | 16.64456 | 63.22278 | 66.57823 | 66.57823 | 83.22278 | 63.22278 | 83.22278 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1954 | 69.26258 | 17.31565 | 66.57823 | 69.26258 | 69.26258 | 86.57823 | 66.57823 | 86.57823 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1955 | 71.41007 | 17.85252 | 69.26258 | 71.41007 | 71.41007 | 89.26258 | 69.26258 | 89.26258 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1956 | 73.12805 | 18.28201 | 71.41007 | 73.12805 | 73.12805 | 91.41007 | 71.41007 | 91.41007 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1957 | 74.50244 | 18.62561 | 73.12805 | 74.50244 | 74.50244 | 93.12805 | 73.12805 | 93.12805 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1958 | 75.60195 | 18.90049 | 74.50244 | 75.60195 | 75.60195 | 94.50244 | 74.50244 | 94.50244 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1959 | 76.48156 | 19.12039 | 75.60195 | 76.48156 | 76.48156 | 95.60195 | 75.60195 | 95.60195 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1960 | 77.18525 | 19.29631 | 76.48156 | 77.18525 | 77.18525 | 96.48156 | 76.48156 | 96.48156 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1961 | 77.74820 | 19.43705 | 77.18525 | 77.74820 | 77.74820 | 97.18525 | 77.18525 | 97.18525 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1962 | 78.19856 | 19.54964 | 77.74820 | 78.19856 | 78.19856 | 97.74820 | 77.74820 | 97.74820 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1963 | 78.55885 | 19.63971 | 78.19856 | 78.55885 | 78.55885 | 98.19856 | 78.19856 | 98.19856 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1964 | 78.84708 | 19.71177 | 78.55885 | 78.84708 | 78.84708 | 98.55885 | 78.55885 | 98.55885 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1965 | 79.07766 | 19.76942 | 78.84708 | 79.07766 | 79.07766 | 98.84708 | 78.84708 | 98.84708 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1966 | 79.26213 | 19.81553 | 79.07766 | 79.26213 | 79.26213 | 99.07766 | 79.07766 | 99.07766 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1967 | 79.40970 | 19.85243 | 79.26213 | 79.40970 | 79.40970 | 99.26213 | 79.26213 | 99.26213 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1968 | 79.52776 | 19.88194 | 79.40970 | 79.52776 | 79.52776 | 99.40970 | 79.40970 | 99.40970 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1969 | 79.62221 | 19.90555 | 79.52776 | 79.62221 | 79.62221 | 99.52776 | 79.52776 | 99.52776 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1970 | 79.69777 | 19.92444 | 79.62221 | 79.69777 | 79.69777 | 99.62221 | 79.62221 | 99.62221 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1971 | 79.75821 | 19.93955 | 79.69777 | 79.75821 | 79.75821 | 99.69777 | 79.69777 | 99.69777 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1972 | 79.80657 | 19.95164 | 79.75821 | 79.80657 | 79.80657 | 99.75821 | 79.75821 | 99.75821 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1973 | 79.84526 | 19.96131 | 79.80657 | 79.84526 | 79.84526 | 99.80657 | 79.80657 | 99.80657 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1974 | 79.87621 | 19.96905 | 79.84526 | 79.87621 | 79.87621 | 99.84526 | 79.84526 | 99.84526 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1975 | 79.90096 | 19.97524 | 79.87621 | 79.90096 | 79.90096 | 99.87621 | 79.87621 | 99.87621 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1976 | 79.92077 | 19.98019 | 79.90096 | 79.92077 | 79.92077 | 99.90096 | 79.90096 | 99.90096 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1977 | 79.93662 | 19.98415 | 79.92077 | 79.93662 | 79.93662 | 99.92077 | 79.92077 | 99.92077 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1978 | 79.94929 | 19.98732 | 79.93662 | 79.94929 | 79.94929 | 99.93662 | 79.93662 | 99.93662 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1979 | 79.95944 | 19.98986 | 79.94929 | 79.95944 | 79.95944 | 99.94929 | 79.94929 | 99.94929 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1980 | 79.96755 | 19.99189 | 79.95944 | 79.96755 | 79.96755 | 99.95944 | 79.95944 | 99.95944 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1981 | 79.97404 | 19.99351 | 79.96755 | 79.97404 | 79.97404 | 99.96755 | 79.96755 | 99.96755 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1982 | 79.97923 | 19.99481 | 79.97404 | 79.97923 | 79.97923 | 99.97404 | 79.97404 | 99.97404 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1983 | 79.98338 | 19.99585 | 79.97923 | 79.98338 | 79.98338 | 99.97923 | 79.97923 | 99.97923 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1984 | 79.98671 | 19.99668 | 79.98338 | 79.98671 | 79.98671 | 99.98338 | 79.98338 | 99.98338 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1985 | 79.98937 | 19.99734 | 79.98671 | 79.98937 | 79.98937 | 99.98671 | 79.98671 | 99.98671 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1986 | 79.99149 | 19.99787 | 79.98937 | 79.99149 | 79.99149 | 99.98937 | 79.98937 | 99.98937 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1987 | 79.99319 | 19.99830 | 79.99149 | 79.99319 | 79.99319 | 99.99149 | 79.99149 | 99.99149 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1988 | 79.99456 | 19.99864 | 79.99319 | 79.99456 | 79.99456 | 99.99319 | 79.99319 | 99.99319 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1989 | 79.99564 | 19.99891 | 79.99456 | 79.99564 | 79.99564 | 99.99456 | 79.99456 | 99.99456 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1990 | 79.99652 | 19.99913 | 79.99564 | 79.99652 | 79.99652 | 99.99564 | 79.99564 | 99.99564 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1991 | 79.99721 | 19.99930 | 79.99652 | 79.99721 | 79.99721 | 99.99652 | 79.99652 | 99.99652 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1992 | 79.99777 | 19.99944 | 79.99721 | 79.99777 | 79.99777 | 99.99721 | 79.99721 | 99.99721 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1993 | 79.99822 | 19.99955 | 79.99777 | 79.99822 | 79.99822 | 99.99777 | 79.99777 | 99.99777 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1994 | 79.99857 | 19.99964 | 79.99822 | 79.99857 | 79.99857 | 99.99822 | 79.99822 | 99.99822 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1995 | 79.99886 | 19.99971 | 79.99857 | 79.99886 | 79.99886 | 99.99857 | 79.99857 | 99.99857 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1996 | 79.99909 | 19.99977 | 79.99886 | 79.99909 | 79.99909 | 99.99886 | 79.99886 | 99.99886 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1997 | 79.99927 | 19.99982 | 79.99909 | 79.99927 | 79.99927 | 99.99909 | 79.99909 | 99.99909 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1998 | 79.99942 | 19.99985 | 79.99927 | 79.99942 | 79.99942 | 99.99927 | 79.99927 | 99.99927 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 1999 | 79.99953 | 19.99988 | 79.99942 | 79.99953 | 79.99953 | 99.99942 | 79.99942 | 99.99942 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2000 | 79.99963 | 19.99991 | 79.99953 | 79.99963 | 79.99963 | 99.99953 | 79.99953 | 99.99953 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2001 | 79.99970 | 19.99993 | 79.99963 | 79.99970 | 79.99970 | 99.99963 | 79.99963 | 99.99963 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2002 | 79.99976 | 19.99994 | 79.99970 | 79.99976 | 79.99976 | 99.99970 | 79.99970 | 99.99970 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2003 | 79.99981 | 19.99995 | 79.99976 | 79.99981 | 79.99981 | 99.99976 | 79.99976 | 99.99976 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2004 | 79.99985 | 19.99996 | 79.99981 | 79.99985 | 79.99985 | 99.99981 | 79.99981 | 99.99981 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2005 | 79.99988 | 19.99997 | 79.99985 | 79.99988 | 79.99988 | 99.99985 | 79.99985 | 99.99985 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2006 | 79.99990 | 19.99998 | 79.99988 | 79.99990 | 79.99990 | 99.99988 | 79.99988 | 99.99988 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2007 | 79.99992 | 19.99998 | 79.99990 | 79.99992 | 79.99992 | 99.99990 | 79.99990 | 99.99990 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2008 | 79.99994 | 19.99998 | 79.99992 | 79.99994 | 79.99994 | 99.99992 | 79.99992 | 99.99992 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2009 | 79.99995 | 19.99999 | 79.99994 | 79.99995 | 79.99995 | 99.99994 | 79.99994 | 99.99994 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| 2010 | 79.99996 | 19.99999 | 79.99995 | 79.99996 | 79.99996 | 99.99995 | 79.99995 | 99.99995 | 0.6 | 0.4 | 0.2 | 20 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
kable(round(t(cbind(
data1$baseline[c(1,2,20,40,66),c("iter block 1")],
data2$baseline[c(1,2,20,40,66),c("iter block 1")],
data3$baseline[c(1,2,20,40,66),c("iter block 1")]
)),digits=3))| 1945 | 1946 | 1964 | 1984 | 2010 |
|---|---|---|---|---|
| 0 | 293 | 218 | 175 | 119 |
| 0 | 89 | 15 | 2 | 1 |
| 0 | 10 | 10 | 10 | 10 |
kable(round(t(dataex$baseline[c(1,2,20,40,66),
c("iter block 1","iter block 2","iter block 3",
"iter block 4","iter block 5","iter block 6",
"iter block 7")]),digits=3))| 1945 | 1946 | 1964 | 1984 | 2010 | |
|---|---|---|---|---|---|
| iter block 1 | 0 | 2 | 2 | 2 | 2 |
| iter block 2 | 0 | 2 | 2 | 2 | 2 |
| iter block 3 | 0 | 2 | 2 | 2 | 2 |
| iter block 4 | 0 | 2 | 2 | 2 | 2 |
| iter block 5 | 0 | 2 | 2 | 2 | 2 |
| iter block 6 | 0 | 2 | 2 | 2 | 2 |
| iter block 7 | 0 | 2 | 2 | 2 | 2 |
Flows in the UK Economy 2014, source: Burgess et al. (2016)
Graph representation of SIM and SIMEX, source: author’s computation
Graph representation of Burgess et al. (2016)
Graph representation of Burgess et al. (2016)
Graph representation of Burgess et al. (2016)
Graph representation of Burgess et al. (2016)
Shiny application
| Households | Production | Government | Central Bank | Sum | |
|---|---|---|---|---|---|
| Money | +H | -H | 0 | ||
| Bills | +Bh | -B | +Bcb | 0 | |
| Net worth | -V | +V | 0 | ||
| Sum | 0 | 0 | 0 | 0 | 0 |
| Households | Production | Government | Central Bank | Sum | |
|---|---|---|---|---|---|
| Consumption | -C | +C | 0 | ||
| Gov. Exp. | +G | -G | 0 | ||
| Income = GDP | +Y | -Y | 0 | ||
| Interests | +r(-1)*Bh(-1) | -r(1)*B(-1) | +r(-1)*Bcb(-1) | 0 | |
| CB profits | +r(-1)*Bcb(-1) | -r(1)*Bcb(-1) | 0 | ||
| Taxes | -T | +T | 0 | ||
| Change in Money | -\(\Delta\) H | +\(\Delta\) H | 0 | ||
| Change in Bills | -\(\Delta\) Bh | +\(\Delta\) B | -\(\Delta\) Bcb | 0 | |
| Sum | 0 | 0 | 0 | 0 | 0 |
#######Determination of output - eq. 4.1
y = cons + g
#######Disposable income - eq. 4.2
yd = y - t + r(-1)*b_h(-1)
#######Tax payments - eq. 4.3
t = theta*(y + r(-1)*b_h(-1))
#######Wealth accumulation - eq. 4.4
v = v(-1) + (yd - cons)
#######Consumption function - eq. 4.5
cons = alpha1*yd + alpha2*v(-1)
#######Cash money - eq. 4.6
h_h = v - b_h
#######Demand for government bills - eq. 4.7
b_h = v*(lambda0 + lambda1*r - lambda2*(yd/v))
#######Supply of government bills - eq. 4.8
b_s = b_s(-1) + (g + r(-1)*b_s(-1)) - (t + r(-1)*b_cb(-1))
#######Supply of cash - eq. 4.9
h_s = h_s(-1) + b_cb - b_cb(-1)
#######Government bills held by the central bank - eq. 4.10
b_cb = b_s - b_h
#######Interest rate as policy instrument - eq. 4.11
r = r_barHow to allocate wealth between bills and money?
Merging the two into Tobin’s portfolio equations:
\[\begin{align} \frac{H_h}{V}&=(1-\lambda_0)-\lambda_1\cdot r+\lambda_2\cdot\left(\frac{YD}{V}\right)\tag{4.6A}\\ \frac{B_h}{V}&=\lambda_0+\lambda_1\cdot r-\lambda_2\cdot\left(\frac{YD}{V}\right)\tag{4.7}\\ H_h&=V-B_h\tag{4.6} \end{align}\]Before doing any modelling, we need to load the package in the R environment.
library(PKSFC)Then, you need to download the attached ‘PC.txt’ file and save it in the folder of your choice. Make sure to set the working directory where you saved the downloaded file. In command line this looks like this but if you use Rstudio, you can use the graphical interface as well (Session>Set Working Directory>Choose Directory)
setwd("pathToYourDirectory")The first thing to do is to load the model and check for completeness.
pc<-sfc.model("Models/PC.txt",modelName="Portfolio Choice Model")
pc<-sfc.check(pc,fill=FALSE)Let’s have a look at the graphical representation of model PC. In order to do so, we need to load the Rgraphviz library:
library("Rgraphviz")We can now look at the graph of model PC:
plot.dag(pc,main="PC" )You can see that there is a cycle in the graph, implying that GDP, taxes, disposable income and consumption are determined all together (and that they fully adapt to any shock applied to the economy). While this is a mathematical property of the system of equations representing the economy we wish to model, it has an economic meaning and you want to be sure that this is what you believe to be the best representation of what you have in mind.
We are now ready to simulate the model
datapc<-simulate(pc)The first of experiment is meant to observe the buffer role of money, in case of random shocks applied to expected disposable income. To do these, we need to modify slightly model pc and change the equations determining consumption, demand for bonds and money, and the expectations on income and wealth. We will call the new model pcRand.
pcRand<-sfc.editEqus(pc,list(list(var="cons",eq="alpha1*yde+alpha2*v(-1)"),
list(var="b_h",eq="ve*(lambda0 + lambda1*r - lambda2*(yde/ve))")))
pcRand<-sfc.addEqus(pcRand,list(
list(var="yde",eq="yd(-1)*(1+rnorm(1,sd=0.1))",
desc="Expected disposable income depending on random shocks"),
list(var="h_d",eq="ve-b_h",desc="Money demand"),
list(var="ve",eq="v(-1)+yde-cons",desc="Expected disposable income")))
pcRand<-sfc.editVar(pcRand,var="yd",init=86.48648)
pcRand<-sfc.check(pcRand)You probably have noticed that the yde equation contains the R function rnorm which allows you to extract a random number from a normal distribution centred around 0, with standard variation of 0.1. The package indeed allows you to use any R function such as min, max or rnorm. However, because of the way the Gauss-Seidel algorithm works, we are facing a problem. Indeed, as the rnorm function extract a number in each iteration, this will mean that the algorithm cannot converge since the difference between each iteration of the algorithm depends of the difference between each extraction. This is why we need to restrict the number of iteration of the Gauss-Seidel.
Before simulating the model, let’s have a look at how the graph of the model has changed:
plot.dag(pcRand,main="PC Random" )You can now see that the cycle observed in the original model has disappeared. Indeed, consumption doesn’t depend on disposable income any more but on expected disposable income. Let’s now simulate the model.
datapcRand<-simulate(pcRand,maxIter=2)This replicates figure 4.1, page 110
plot(pcRand$time,datapcRand$baseline[,"h_h"],type="l",xlab="",ylab="",lty=1,lwd=2,
ylim=c(1*min(datapcRand$baseline[,c("h_h","h_d")],na.rm=T),
1.2*max(datapcRand$baseline[,c("h_h","h_d")],na.rm=T)))
lines(pcRand$time,datapcRand$baseline[,"h_d"],lty=2,lwd=2)
legend(x=1950,y=1.2*max(datapcRand$baseline[,c("h_h","h_d")],na.rm=T),
legend=c("Money held","Money demand"),lty=c(2,1),lwd=2,bty="n")This graph highlights the buffer role of certain stocks in PK-SFC models. Indeed, because expectations are incorrect or because demand might not be equal to supply in any market, at least one stock will not be equal to the targeted level. As highlighted by Foley (1975), in a model without perfect foresight you need a buffer stock in order to obtain equilibrium between demand and supply. The role of buffer stocks in PK-SFC model is thus fundamental and is at the hart of the approach used by Godley (1999) in his seven unsustainable processes. It is by observing the dynamics of certain stock-flow norms that you are able to observe the unsustainable processes evolving in an economy, because stocks precisely absorb disequilibrium.
The graphs presented on the paper are based on the model PCEX, which includes expectation on disposable income and wealth. We will first create the model before simulating it.The shock represents an 100 basis point increase in interest rates.
pcex<-sfc.editEqus(pc,list(
list(var="cons",eq="alpha1*yde+alpha2*v(-1)"),
list(var="b_h",eq="ve*(lambda0 + lambda1*r - lambda2*(yde/ve))")))
pcex<-sfc.addEqus(pcex,list(
list(var="yde",eq="yd(-1)",desc="Expected disposable income"),
list(var="h_d",eq="ve-b_h",desc="Money demand"),
list(var="ve",eq="v(-1)+yde-cons",desc="Expected disposable income")))
pcex<-sfc.editVar(pcex,var="yd",init=86.48648)
pcex<-sfc.check(pcex)
init = datapc$baseline[56,]
pcex<-sfc.addScenario(model=pcex,vars="r_bar",values=0.035,inits=1960,
ends=2010,starts=init)
datapcex<-simulate(pcex)This replicates plot figure 4.3. p. 112
time2=c(1950:2000)
plot(time2, datapcex$scenario[as.character(time2),"h_h"]/
datapcex$scenario[as.character(time2),"v"],
type="l",xlab="",ylab="",lty=1,lwd=2)
par(new=T)
plot(time2,datapcex$scenario[as.character(time2),"b_h"]/
datapcex$scenario[as.character(time2),"v"],
lty=2, lwd=2,type="l",axes=F,ylab="",xlab="")
axis(4,pretty(c(0.750, 1.1*max(datapcex$scenario[,"b_h"]/
datapcex$scenario[,"v"],na.rm=T))))
legend(x=1970,y=0.78,legend=c("Share of Bills",
"Share of Money"),lty=c(2,1),lwd=2,bty="n")
grid()No surprise here, the change in return rates coming from the shock led to a re-allocation between Bills and Money, in favour of the former.
This replicates fig 4.4 . p 113
time2=c(1950:2000)
matplot(time2,datapcex$scenario[as.character(time2),c("yd","cons")]
,type="l",xlab="",ylab="",lty=1:2,col=1,lwd=2)
legend(x=1970,y=88,legend=c("Disposable Income","Consumption"),
lwd=c(2,2),lty=c(1,2),bty="n")
grid()These results are more surprising as increased interest rates lead to an increase both in disposable income and consumption. To understand this, one has to bear in mind that increase in interest rates leads to an increase in government spending and thus to an increase in the fiscal stance which determines the steady state:
\(Y^\star=\frac{G+r \cdot B^{\star}_h \cdot (1-\theta)}{\theta}\)
We change the propensities to consume to incorporate the fact that they might be impacted by interest rates and thus have
\[\begin{equation} \alpha_1=\alpha_10-\iota\cdot r_{-1} \tag{4.32} \end{equation}\]This piece of code shows how to implement this
pc_end<-sfc.addEqus(pcex,list(
list(var="alpha1",eq="alpha10-iota*r(-1)")))
pc_end<-sfc.addVars(pc_end,list(
list(var="alpha10",init="0.7",
desc="Endoegnous propensity to consume - autonomous term"),
list(var="iota",init="10",
desc="Endoegnous propensity to consume - interest rate impact")))
pc_end$scenarios<-NULL
pc_end<-sfc.check(pc_end)
datapc_end<-simulate(pc_end)
init = datapc_end$baseline[56,]
pc_end<-sfc.addScenario(model=pc_end,vars="r_bar",values=0.035,
inits=1960,ends=2010,starts=init)
datapc_end<-simulate(pc_end)This replicates plot figure 4.9. p. 123
time2=c(1950:2000)
matplot(time2,
datapc_end$scenario_1[as.character(time2),c("v","yd","cons","y")],
type="l",xlab="",ylab="",lty=1:4,lwd=2,col=1)
grid()
legend(x=1950,y=100,legend=c("Households wealth","Disposable income",
"Consumption","GDP"),lty=1:4,lwd=2,bty="n")This shows that now the short-term dynamics display the Keynesian paradox of thrift but the long-run steady state still shows a positive impact of interests rates on GDP. The short-run recession is governed by the decrease in the propensity to consumed and is slowly compensated by the increase in wealth (and the increased public debt as a counterpart) leading to the new steady state. The steady state is still driven by the fiscal stance.
This replicates plot figure 4.10. p. 123
timelag=c(1950:2000)
time2=c(1951:2001)
plot(time2,
datapc_end$scenario_1[as.character(time2),c("g")]+
datapc_end$scenario_1[as.character(timelag),c("r")]*
datapc_end$scenario_1[as.character(timelag),c("b_s")]-
datapc_end$scenario_1[as.character(timelag),c("r")]*
datapc_end$scenario_1[as.character(timelag),c("b_cb")],
type="l",xlab="",ylab="",lty=1:4,lwd=2,col=1,ylim=c(20.8,23))
lines(time2,datapc_end$scenario_1[as.character(time2),"t"],type="l",xlab="",
ylab="",lty=2,lwd=2)
grid()
legend("topleft",legend=c("Government expenditures plus net debt services",
"Tax revenues"),lty=1:4,lwd=2,bty="n")The Maastricht treaty: The reference values referred to […] are:
The value of a financial asset is supposed to reflect the net present value of future cash flows. Thus for a long-term bond who is never redeemed, the price \(p_{bL}\) of the bond is given by
\(p_{bL} = \sum \frac{1}{(1+r_{bL})^t} = \frac{1}{r_{bL}}\)
Return rate of an asset is equal to both the yield (interests, dividends) and the capital gains normalized to the nominal value of the asset. In the case of a long-term bond:
\(Rr_{bL} = r_{bL}(-1)+\frac{\Delta p_{bL}}{p_{bL}(-1)} = \frac{1+p_{bL}-p_{bL}(-1)}{p_{bL}(-1)}\)
However, what matter for the households when making their decision is the expected price of bonds \(p_{bL}^e\), given the current price of bonds. This leads to the pure expected rate of return.
\(PERr_{bL} = r_{bL} + \frac{p^e_{bL}-p_{bL}}{p_{bL}}\)
Because expectation are not followed by everyone, households might incorporate a weight \(\xi\) into their expectation formation. This leads to the expected return rate
\(ERr_{bL} = r_{bL} + \xi \frac{p^e_{bL}-p_{bL}}{p_{bL}}\)
The change in prices of bonds leads to capital gains which is one of two parts of the change in nominal value of financial assets. More precisely:
\[\begin{align*} \Delta (p_{BL}\cdot BL) &=\left(p_{bL} \cdot BL\right) - \left(p_{bL}(-1) \cdot BL(-1)\right)\\ &= p_{BL} \cdot BL-p_{BL}\cdot BL(-1) + p_{BL}\cdot BL(-1)-p_{bL}(-1) \cdot BL(-1)\\ &= p_{BL} \left( BL-BL(-1)\right) + BL(-1) \left(p_{BL}-p_{BL}(-1)\right)\\ &= p_{BL} \cdot \Delta BL+ BL(-1)\cdot \Delta p_{BL} \end{align*}\]Thus total changes in nominal values of an asset is equal to the quantity of assets bought (or sold) at the current price (\(p_{BL} \cdot \Delta BL\)) and to the capital gains (or losses) made on the stock of assets held at the beginning of the period (\(BL(-1)\cdot \Delta p_{BL}\))
Ostergaard diagram
| Households | Production | Government | Central Bank | Sum | |
|---|---|---|---|---|---|
| Money | +H | -H | 0 | ||
| Bills | +Bh | -B | +Bcb | 0 | |
| Bonds | +BL.pbL | -BL.pbL | 0 | ||
| Net worth | -V | +V | 0 | ||
| Sum | 0 | 0 | 0 | 0 | 0 |
| Households | Production | Government | Central Bank | Sum | |
|---|---|---|---|---|---|
| Consumption | -C | +C | 0 | ||
| Gov. Exp. | +G | -G | 0 | ||
| Income = GDP | +Y | -Y | 0 | ||
| Interests on bills | +rb(-1)*Bh(-1) | -rb(1)*B(-1) | +rb(-1)*Bcb(-1) | 0 | |
| Interests on bonds | +BL(-1) | -BL(-1) | 0 | ||
| CB profits | +r(-1)*Bcb(-1) | -r(1)*Bcb(-1) | 0 | ||
| Taxes | -T | +T | 0 | ||
| Change in Money | -\(\Delta\) H | +\(\Delta\) H | 0 | ||
| Change in Bills | -\(\Delta\) Bh | +\(\Delta\) B | -\(\Delta\) Bcb | 0 | |
| Change in Bonds | -\(\Delta\) BL.pbL | +\(\Delta\) BL.pbL | 0 | ||
| Sum | 0 | 0 | 0 | 0 | 0 |
| Memo: Capital gains | -\(\Delta\) pbL.BL(-1) | +\(\Delta\) pbL.BL(-1) | 0 |
# MODEL
# Determination of output - eq. 5.1
y = cons + g
# Regular disposable income - eq. 5.2
yd_r = y - t + r_b(-1)*b_h(-1) + bl_h(-1)
# Tax payments - eq. 5.3
t = theta*(y + r_b(-1)*b_h(-1) + bl_h(-1))
# Wealth accumulation - eq. 5.4
v = v(-1) + (yd_r - cons) + cg
# Capital gains on bonds - eq. 5.5
cg = (p_bl-p_bl(-1))*bl_h(-1)
# Consumption function - eq. 5.6
cons = alpha1*yd_r_e + alpha2*v(-1)
# Expected wealth - eq. 5.7
v_e = v(-1) + (yd_r_e - cons) + cg
# Cash money - eq. 5.8
h_h = v - b_h - p_bl*bl_h
# Demand for cash - eq. 5.9
h_d = v_e - b_d - p_bl*bl_d
# Demand for government bills - eq. 5.10
b_d = v_e*(lambda20 + lambda22*r_b - lambda23*er_rbl - lambda24*(yd_r_e/v_e))
# Demand for government bonds - eq. 5.11
bl_d = v_e*(lambda30 - lambda32*r_b + lambda33*er_rbl - lambda34*(yd_r_e/v_e))/p_bl
# Bills held by households - eq. 5.12
b_h = b_d
# Bonds held by households - eq. 5.13
bl_h = bl_d
# Supply of government bills - eq. 5.14
b_s = b_s(-1) + (g + r_b(-1)*b_s(-1) + bl_s(-1)) - (t + r_b(-1)*b_cb(-1)) -p_bl*(bl_s-bl_s(-1))
# Supply of cash - eq. 5.15
h_s = h_s(-1) + b_cb - b_cb(-1)
# Government bills held by the central bank - eq. 5.16
b_cb = b_s - b_h
# Supply of government bonds - eq. 5.17
bl_s = bl_h
# Expected rate of return on bonds - eq. 5.18
er_rbl = r_bl+chi*(p_bl_e - p_bl)/p_bl
# Interest rate on bonds - eq. 5.19
r_bl = 1/p_bl
# Expected price of bonds - eq. 5.20
p_bl_e = p_bl
# Expected capital gains - eq. 5.21
cg_e = chi*(p_bl_e - p_bl)*bl_h
# Expected regular disposable income - eq. 5.22
yd_r_e = yd_r(-1)
# Interest rate on bills - eq. 5.23
r_b = r_bar
# Price of bonds - eq. 5.24
p_bl = p_bl_barLet’s have a look at the graph of model LP
lp<-sfc.model(fileName="Models/LP.txt",modelName="Model Liquidity Preference")
plot.dag(lp)The vertical conditions state that the sum of the exogenous components sums to 1, i.e. you cannot allocate more than the wealth and that the effect of the various return rates and disposable income are subtracting the holding of one asset to allow the increase of the holding of another.
\[\begin{align} \lambda_{10} + \lambda_{20} + \lambda_{30} &= 1 \tag{ADUP.1}\\ \lambda_{11} + \lambda_{21} + \lambda_{31} &= 0 \tag{ADUP.2}\\ \lambda_{12} + \lambda_{22} + \lambda_{32} &= 0 \tag{ADUP.3}\\ \lambda_{13} + \lambda_{23} + \lambda_{33} &= 0 \tag{ADUP.4}\\ \lambda_{14} + \lambda_{24} + \lambda_{34} &= 0 \tag{ADUP.5} \end{align}\]The horizontal conditions, added by Godley (1996) state that the impact of an increase in the ‘own rate’ of an asset should be equal to the impact of a fall in all the other rates.
\[\begin{align} \lambda_{11} &= - \left( \lambda_{12} + \lambda_{13}\right) \tag{ADUP.6}\\ \lambda_{22} &= - \left( \lambda_{21} + \lambda_{23}\right) \tag{ADUP.7}\\ \lambda_{33} &= - \left( \lambda_{31} + \lambda_{32}\right) \tag{ADUP.8} \end{align}\]Other authors propose the symmetry conditions that state that an increase in the return rate of one asset A will have the same impact on the holding of asset B as the increase in the return rate of asset B will have on the holding of asset A
\[\begin{align} \lambda_{12} &= \lambda_{21} \tag{ADUP.9}\\ \lambda_{13} &= \lambda_{31} \tag{ADUP.10}\\ \lambda_{23} &= \lambda_{32} \tag{ADUP.11} \end{align}\]Note that vertical + symmetry imply automatically horizontal but that vertical + horizontal does not necessarily imply symmetry.
The first scenario run looks at the impact of interest rates (short- and long-term) on real demand. We assume that the government increases short-term rates from 3% to 4% and the long-term rate from 5% to 6.66%. We assume that this increase is completely unexpected and that the Treasury is believed in that no other change will occur.
lp<-sfc.addScenario(model=lp,vars=list(c("r_bar","p_bl_bar")),values=list(c(0.04,15)),
inits=1960,ends=2000)
datalp<-simulate(lp)This replicates plot figure 5.2. p. 152
time2=c(1950:2000)
plot(time2, datalp$scenario_1[as.character(time2),"v"]/
datalp$scenario[as.character(time2),"yd_r"],
type="l",xlab="",ylab="",lty=1,lwd=2)
legend("bottomright",legend=c("Wealth to income ratio"),lty=c(1),lwd=2,bty="n")
grid()We see that the increase in long-term interest rate leads to capital losses on bonds, leading to a decrease in wealth to income ratio. However because interest rates have increase, households have now a higher disposable income leading to a replenishing of their wealth. Note that as in model PC, the steady state GDP is given by
\(Y^\star=\frac{G+\left(r_B\cdot B^{\star}_h+BL_S^\star\right) \cdot (1-\theta)}{\theta} = \frac{G_{NT}}{\theta}\)
This replicates plot figure 5.3. p. 152
time2=c(1950:2000)
matplot(time2, datalp$scenario_1[as.character(time2),c("yd_r","cons")],
type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("topleft",legend=c("Disposable income","Consumption"),lty=c(1,2),lwd=2,bty="n")
grid()This replicates plot figure 5.4. p. 153
time2=c(1950:2000)
matplot(time2, cbind(datalp$scenario_1[as.character(time2),c("b_h")],
datalp$scenario_1[as.character(time2),c("p_bl")]*
datalp$scenario_1[as.character(time2),c("bl_h")])/
datalp$scenario_1[as.character(time2),c("v")],
type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("center",legend=c("Bonds to wealth ratio","Bills to wealth ratio"),
lty=c(2,1),lwd=2,bty="n")
grid()In order to introduce liquidity preferences, we need to modify the model. The first modification is to change the expectation on bonds price to include (i) a learning element and (ii) a stochastic term. Because expectation are going to play a role, we also need to modify the price equation for bonds such that it now reflects a corridor approach to bond pricing. The assumption is thus that the Treasury will try to pin down the price of bonds but only to a certain extent. If the demand is too large (or too low) the treasury will let prices (and interest) float. We assume that the government aims at a certain debt structure in terms of maturity.
#Adding the equations
lp2<-sfc.addEqus(lp,list(
list(var="z1",equ = "TP>top"),
list(var="z2",equ = "TP<bot"),
list(var="TP",equ = "(bl_h(-1)*p_bl(-1))/(bl_h(-1)*p_bl(-1)+b_h(-1))")))
#Adding the parameters
lp2<-sfc.addVars(lp2,list(
list(var="betae",init=0.5),
list(var="beta",init=0.01),
list(var="add",init=0),
list(var="top",init=0.495),
list(var="bot",init=0.505)
))
#Adding initial values to the endogenous variables that have lags
lp2<-sfc.editEnd(lp2,var="p_bl_e", init=20)
lp2<-sfc.editEnd(lp2,var="p_bl", init=20)
#Modifying the existing equations
lp2<-sfc.editEqus(lp2,list(
list(var="p_bl", eq="(1+z1*beta-z2*beta)*p_bl(-1)"),
list(var="p_bl_e",equ = "p_bl_e(-1)-betae*(p_bl_e(-1) - p_bl) + add")))
#Removing all existing scenarios (coming from LP)
lp2$scenarios<-NULL
#Checking that the model is complete
lp2<-sfc.check(lp2,fill=F)
#Adding the scenario
lp2<-sfc.addScenario(model=lp2,vars=list(c("r_bar")),values=list(c(0.035)),inits=1960,
ends=2000)
lp2<-sfc.addScenario(model=lp2,vars=list(c("add")),values=list(c(-3)),inits=c(1955),
ends=1955)Let’s have a look at the graph of model LP2
plot.dag(lp2)We are now ready to simulate the model
datalp2<-simulate(lp2)This replicates plot figure 5.5. p. 156
time2=c(1950:2000)
plot(time2,datalp2$scenario_1[as.character(time2),"r_b"],
lty=2, lwd=2,type="l",ylab="",xlab="",ylim=c(0.029,0.036))
par(new=T)
plot(time2, datalp2$scenario_1[as.character(time2),"r_bl"],
type="l",xlab="",axes=F,ylab="",lty=1,lwd=2,ylim=c(0.0495,0.052))
axis(4,pretty(c(0.0495,0.052)))
legend("bottomright",legend=c("Long term interest rate",
"Short term interest rate"),lty=c(1,2),lwd=2,bty="n")
grid()This replicates plot figure 5.6. p. 156
time2=c(1950:2000)
matplot(time2, datalp2$scenario_1[as.character(time2),c("TP","bot","top")],
type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend(x=1950,y=0.502,legend=c("Share of bonds in government debt held by households"),
lty=c(2),lwd=2,bty="n")
grid()Let’s add a second scenario by which there is an expected fall in the price of long-term bonds.
lp2<-sfc.addScenario(model=lp2,vars=list(c("add")),values=list(c(-3)),inits=c(1955),
ends=1955)
datalp2<-simulate(lp2)This replicates plot figure 5.7. p. 157
time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("r_bl")],
type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend(x=1950,y=0.502,legend=c("Share of bonds in government debt held by households"),
lty=c(2),lwd=2,bty="n")
grid()This replicates plot figure 5.8. p. 158
time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("p_bl_e","p_bl")],
type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("center",legend=c("Expected price of bonds","Actual price of bonds"),
lty=c(1,2),lwd=2,bty="n")
grid()This replicates plot figure 5.9. p. 158
time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("TP","bot","top")],
type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend("center",legend=c("Share of bonds in government debt held by households"),
lty=c(2),lwd=2,bty="n")
grid()#Supply of consumption goods - eq. 7.1
c_s = c_d
#Supply of investment goods - eq. 7.2
i_s = i_d
#Supply of labour - eq. 7.3
n_s = n_d
#Transactions of the firms
#GDP - eq. 7.5
y = c_s + i_s
#Wage bill - eq. 7.6
wb_d = y - r_l(-1)*l_d(-1) - af
#Depreciation allowances - eq. 7.7
af = delta*k(-1)
#Demand for bank loans - eq. 7.8
l_d = l_d(-1) + i_d - af
#Transactions of households
#Disposable income - eq. 7.9
yd = wb_s + r_m(-1)*m_h(-1)
#Bank deposits held by households - eq. 7.10
m_h = m_h(-1) + yd - c_d
#The wage bill
#"Supply" of wages - eq. 7.13
wb_s = w*n_s
#Labour demand - eq. 7.14
n_d = y/pr
#Wage rate - eq. 7.15
w = wb_d/n_d
#Household behaviour
#Demand for consumption goods - eq. 7.16
c_d = alpha0 + alpha1*yd + alpha2*m_h(-1)
#The investment behaviour
#Accumulation of capital - eq. 7.17
k = k(-1) + i_d - da
#Depreciation allowances - eq. 7.18
da = delta*k(-1)
#Capital stock target - eq. 7.19
k_t = kappa*y(-1)
# below is an additional equation I've defined to get output to capital ratio as in figure 7.4, the model will run even without equation
OCR = y/k(-1)
#Demand for investment goods - eq. 7.20
i_d = gamma*(k_t - k(-1)) + daBMW<-sfc.model("Models/BMW.txt")## Warning in sfc.check(model, fill = fill): The following variables have lags
## but no initial values: - r_l - r_m
plot.dag(BMW)| Components | Firms Current account | Firms Capital account |
|---|---|---|
| Sales | \(+S\) | |
| Change in the value of inventories | +\(\Delta IN\) | -\(\Delta IN\) |
| Wages | \(-WB\) | |
| Interest on loans | \(-r_{l-1} L_{-1}\) | |
| Entrepreneurial profits | \(-F\) | |
| Change in loans | \(+\Delta L\) | |
| Sum | 0 | 0 |
Foley, D. 1975. “On Two Specifications of Asset Equilibrium in Macroeconomic Models.” Journal of Political Economy 83-2 (2): 303–24.
Godley, Wynne. 1999. “Seven Unsustainable Processes: Medium-Term Prospects and Policies for the United States and the World.” Strategic Analysis. The Levy Economic Institute of Bard College.